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ABSTRACT 

Optical properties of hybrid plasmonic waveguides and of low-Q cavities, formed by waveguides of finite length 
are investigated numerically. These structures are of interest as building-blocks of plasmon lasers. We use 
a time- harmonic finite-element package including a propagation- mode solver, a resonance- mode solver and a 
scattering solver for studying various properties of the system. Numerical convergence of all used methods is 
demonstrated. 
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Figure 1. Visualization of a plasmonic waveguide mode at wavelength Ao = 489 nm. Left: Projection of the vectorial 
electric field onto the x^/-plane. Right: Logarithm of the magnitude of the y component of the electric field (black/dark: 
low intensity, yellow/bright: high intensity). The total field of view has a size of 90nmx95nm. 



1. INTRODUCTION 

Recently, plasmon lasers, sometimes also called nanolasers or spasers (surface plasmon lasers) P have been realized 
in different systems. Essentially, in these devices light can be confined at deep subwavelength scale. This opens 
new prospects for photonic integration and applications like biosensing on a nanoscaleP Approaches for light 
confinement include nanosphere arrangements P rectangular semiconductor nanopillars encapsulated in silver P 
semiconductor nanowires on top of a silver substrate,^ bowtie-shaped silver nanoparticle cavities,^ and others. In 
some of the devices, plasmon-based waveguiding structures play a crucial role. Plasmon waveguides incorporating 
metal and several different dielectric materials are able to confine the light field in two dimensions to a very small 
spot size (see typical mode field distribution in Figure [T]) at relatively low propagation lossP"^ To distinguish 
these waveguides from simpler approaches consisting of only metal and a single dielectric material they are also 
called hybrid plasmonic waveguides P 

Fast and accurate 3D Maxwell solvers are needed for designing structural parameters of such plasmonic 
waveguides. Due to the multi-scale nature of the corresponding field distributions, accurate computation of 
the properties of such devices can be numerically challenging. We have developed finite-element method (FEM) 
based solvers for the Maxwell eigenvalue and for the Maxwell scattering problems. The method is based on higher 
order vectorial elements, adaptive unstructured grids, and on a rigorous treatment of transparent boundaries. 
The method has been applied to plasmonic devices like plasmonic antennas, gratings and waveguides PES Here 
we investigate plasmonic waveguides and cavities in layouts similar to the devices used in experiments by the 
group of X. Zhang. 5 




Figure 2. Left: Schematic of the setup: CdS nanowire of diameter D and length L placed on a thin layer of MgF2 
(hi = 5nm) on a silver layer (subspace). Right: Mesh of the geometry of a plasmonic waveguide with finite length 
L = 374 nm in a 3D setup. 

2. INVESTIGATED SETUP 

The investigated setup consists of a cadmium sulphide (CdS) semiconductor nanowire placed on a thin layer of 
(insulating) magnesium fluoride (MgF2) on a silver (Ag) surface. Figure [2] (left) shows the schematic of such a 
device. The cadmium sulphide nanowire with diameter D extends in z-direction over a length L. The parameters 
chosen in this numerical study are D = 50 nm, hi = 5nm, L of up to few jam. The CdS material is pumped 
optically and lases at a wavelength of A = 489 nm. However, due to the small size of the structure the supported 
propagating modes are concentrated mainly in the MgF2 region below the nanowire, with some overlap to the 
amplifying CdS material. In the x?/-cross-section these modes are concentrated at deep subwavelength scale 
due to the plasmonic nature of the excitation at the Ag - MgF2 interface. The chosen geometrical parameters 
correspond to experimentally realized designs P In our numerical study of optical properties of the plasmonic 
waveguide design we assume the following relative permittivities: metal layer: e r — —9.008 + 0.3i, insulating 
layer: e r = 1.9047, nanowire: e r = 6.3, air: e r — 1. Please note that here we study only the passive optical 
properties of waveguide and cavity (assumption of a real permittivity of the nanowire). 




Figure 3. Left: Plasmonic waveguide geometry (x?/-cross-section through the geometry depicted in Fig. [2}: CdS nanowire 
(blue), MgF2 layer (dark grey), Ag (red), surrounded by air (light grey). Center: Mesh of the geometry for the FEM 
discretization. Right: Adaptively refined mesh (two refinement steps). 



The 3D geometry and corresponding mesh of a plasmonic waveguide with finite length L are shown in Figure [2j 
The geometry is discretized using prismatoidal elements. The subgridding of the 3D mesh in z-direction is not 
shown here. For light scattering simulations, a source field which typically is a 2D waveguide mode can be 
coupled into the computational domain at the domain boundaries. For eigenvalue computations modelling the 
resonance of the cavity, no source fields are applied. Please note that due to a mirror symmetry of the geometry, 
the computational domain can be reduced to the half space x > 0. Figure [3] shows the geometry of the plasmonic 
waveguide x?/-cross-section and the finite-element mesh discretizing the geometry. Here, the geometry is assumed 
to be invariant in z-direction; MgF2, Ag and air regions are modeled as extending to infinity in the xy-plane. 
The computational domain on which the mode field is computed is discretized with triangular elements while 
the surrounding exterior domain on which purely outgoing fields are assumed is discretized with quadrilateral 
elements. The more complex subgridding of the exterior domain for the adaptive realization of PMLs is not 
shown here.^ 

3. NUMERICAL RESULTS 
3.1 Guided modes in a plasmonic waveguide 

The geometry of the plasmonic waveguide is invariant z-direction, cf. Figure [2j A propagating mode is a 
solution to the time harmonic Maxwell's equations with frequency cj, which exhibits a harmonic dependency in 
^-direction: 

E = E pm (>, y) exp (ik z z) . 

E pm (^j2/) is the electric propagation mode on the 2D cross section and the parameter k z is called propagation 
constant. The effective refractive index n e fr is defined as n e ^ = k z /ko with fco = 27r/Ao, where Ao is the 
vacuum wavelength of light. For finding propagating modes we solve an eigenvalue problem for the propagation 
constant k z and propagation mode E pm (x, ?/), using the propagation mode solver included in our programme 
package JCMsuite. 13 Figure [T] shows a computed field distribution for the physical parameters as defined in 
Section [2j The mode is concentrated at a very small cross-section. 7 The computed effective refractive index 
for the given parameters is n e ^ ~ 1.2374 + 0.0118i. We have checked that we can compute this number to a 
high numerical accuracy with relative errors of the real and imaginary parts of below 10 _T . Figure |4] shows the 
corresponding convergence plot: The relative error of n e ff is computed for different levels of grid refinement. The 
higher the level of grid refinement, the higher is the number of unknowns in the FEM problem. The relative error 
of a measured quantity a(N) (here real and imaginary parts of n e ff) is defined as Aa(N) = \ot(N) — ce qe |/|^qe| 7 
where a qe is the quantity computed from the quasi-exact solution, i.e., from a solution on a finer mesh than the 
meshes of the solutions corresponding to a(N). It would be desirable to compare a(N) to an analytical solution, 




Figure 4. Computation of propagating modes in a plasmonic waveguide. Left: Convergence of the relative error of the 
real and imaginary parts of the eigenvalue (propagation constant) with number of unknowns of the FEM problem, N. 
Right: Computational effort in terms of cpu-time with number of unknowns for the same data set. 



but for problems where an analytical solution is not available, the quasi-exact solution is used as a makeshift. In 
previous works we also have checked that our methods converge to analytical results^ and - within the various 
numerical uncertainties - to the same numerical results as other software implementations of FEM and of other 
rigorous methods for solving Maxwell's equationsP^ Here, we have used finite elements of second polynomial 
order; typical computation times on a standard PC are in the range of seconds to minutes for the plotted results. 

3.2 3D simulation of light propagation over a finite distance of an infinite plasmonic 
waveguide 

In this Section results on 3D light scattering simulations of mode propagation through the investigated waveguide 
are presented. We performed simulations of light propagation over a certain distance L of the infinite waveguide. 
The eigenvalue computed with the propagating mode solver, as described in Section [3TTJ provides a very accurate 
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Figure 5. Field distribution of a waveguide mode propagating through the 3D plasmonic waveguide. Top: xz-cross- 
section, ^/-component of the electric field in a pseudocolor visualization. Bottom: 2/z-cross-section, ^/-component of the 
electric field in a pseudocolor visualization. 



value of expected transmission T over a certain distance. The transmission is given by T = exp(— 47rL^(n e ff) / Ao). 
We can use this number to check the accuracy of numerical results for specific numerical parameter settings for 
the 3D setup. 
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Figure 6. Computation of light propagating through a 3D plasmonic waveguide. Convergence of the relative error of 
transmission, AT, over a distance of L — 1.5 /mi with number of unknowns of the FEM problem, N. 

For simulating transmission through the waveguide at specific wavelength Ao, we first compute the funda- 
mental propagation mode of the waveguide at Ao as described in Section [3TTT The obtained mode field is applied 
as input data to one of the boundaries of the 3D computational domain (front boundary in Fig. [2]), such that 
the mode propagates in -f z-direction. We then compute the scattered light field in the setup corresponding 
to this excitation using higher-order finite-elements. A typical 3D computational domain is depicted in Fig. [2j 
Transparent boundary conditions take into account the specific geometry of the problem where waveguides are 
modelled to extend to infinity in the exterior domain.^ In post-processes we extract energy fluxes through 
interfaces and field distributions in several cross-sections from the 3D light field distribution. The transmission 
is given by the quotient of transmitted and incoming energy flux, T = E t /E[ n . 

Figure [5] visualizes the computed electric field intensity in a cross-section in a xz-plane in the center of the 
MgF2 layer (top) and in a cross-section in a ?/z-plane at x = (through the center of the nanowire, bottom). 
Please note the slight decay of the amplitude by a factor of about 0.8 over the propagation distance of L = 1.5 /xm. 
Figure [6] shows how the transmission through the waveguide computed with 3D light scattering simulations 
converges to the quasi-exact result obtained from the 2D propagation mode computation. Accuracies in the 
range of 0.1% are easily reached. Computation times for the results shown in this plot range between one and 
few minutes on a standard workstation with RAM below 10GB. 

3.3 3D simulation of transmission through a finite plasmonic waveguide 

For investigating resonances of a finite plasmonic waveguide we performed a numerical experiment where we 
simulate transmission through such a device. Figure [7] shows a schematic of the setup, a corresponding FEM 
mesh is shown in Figure [8] (left). Coupling to the finite waveguide and outcoupling from the waveguide is done 
by identical, semi-infinite plasmonic waveguides, which are separated from the finite waveguide by a distance 
W. The total length of the computational domain in these simulations is fixed to 2.2 jiva. Figure [71 (right) 
shows transmission spectra through the plasmonic waveguide resonators. In these spectra the resonator length 
L has been varied at fixed vacuum wavelength of the incoming waveguide mode of 489.0 nm. We have recorded 
spectra for various widths W of the gaps separating the resonator from the in- and outcoupling waveguides. 
Clearly, resonance peaks are observed in the transmission spectra. As expected, for increasing gap width W, 
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Figure 7. Left: Schematic of the setup (top view and sideview): CdS nanowire of diameter D and length L placed on a thin 
layer of MgF2 (hi = 5nm) on a silver layer (subspace). Coupling to the nanowire by identical, semi-infinite waveguides 
on either side of the device, separated by a distance W. Right: Transmission T through plasmonic waveguides of varied 
length L, separated from incoupler and outcoupler waveguides by distances W = 100, 150, 200 nm. Vacuum wavelength 
of the incoming waveguide modes is 489.0 nm. 



the overall transmission decreases. Also resonance widths and positions slightly shift with the mirror width 
W. Figure [8] shows field distributions in cross-sections through the 3D computational domain for W = 200 nm 
and L = 370 nm, and for W = 200 nm and L — 1360 nm, both at a transmission maximum. High intensities 
at the end facets of the resonator can be observed (cf. a comparable experimental image by Oulton et alW), 
corresponding to a standing- wave excitation of the resonator. 

3.4 Direct simulation of resonances of a finite plasmonic waveguide 

As has been experimentally demonstrated,^ a plasmonic waveguide of finite length can be used as nanolaser cavity. 
Computing Q-factor and resonance wavelength of such a resonator from transmission spectra of a resonator, as 
obtained from light-scattering simulations (see Section 13. 3p can be time-consuming, and - more importantly - 
the scattering response of a low-Q resonator makes quantitative investigations of the optical properties difficult 
due to the presence of a strong diffuse scattered field. For computing optical properties of such resonator we 
therefore use a FEM eigenmode solver included in JCMsuite. We have recently shown that this solver allows to 
compute resonance wavelength and Q-factor for high-Q photonic crystal microcavities at very high accuracy.^ 




Figure 8. Top left: Mesh of the 3D setup. Top right: Field distribution of a resonance mode for L = 370nra and 
W — 200 nm: xz-and ^/^-cross-sections, electric field intensity in a pseudocolor visualization (logarithmic scale). Bottom: 
L = 1360nra, W = 200 nm, xz-cross-section, electric field intensity in a pseudocolor visualization (logarithmic scale). 
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Figure 9. Left: Dependence of resonance wavelength Ao on cavity length L. Right: Dependence of cavity Q-factor on L. 



Given a geometrical setup with specific material parameters, one computes an electric field distribution E and 
a complex eigenfrequency uj which satisfy Maxwell's time-harmonic wave equation 

V x /i^V x E = uo 2 eE 

without sources; electric permittivity and magnetic permeability are denoted by e and /i, respectively. When 
the eigenmode (E, uj) is computed, the respective Q-factor is deduced from the real and imaginary parts of the 
complex eigenfrequency, Q — 5ft(u;)/(— 29(o;)), the resonance wavelength Ao is given by Ao = 27rco/3?(o;), with 
the speed of light cq. 

In order to find the cavity length L for which the cavity is resonant at the target wavelength Ao = 489 nm, we 
have performed simulations with varied L and otherwise fixed parameters. The resulting plots on the behavior 
of resonance wavelength and Q-factor change with cavity length are shown in Figure [9] For these results, we 
have again chosen physical parameters as defined in Section [2j As can be seen from the Figure, a resonance 




Figure 10. Left: Field distribution of a resonance mode for L — 374nra and Ao ~ 489 nm. Top: xz-cross-section, y- 
component of the electric field in a pseudocolor visualization (logarithmic scale). Bottom: ^z-cross-section, ^-component 
of the electric field in a pseudocolor visualization (logarithmic scale). Right: Numerical convergence of the simulated 
eigenvalues: Dependence of the relative errors of the real and imaginary part of the complex eigenfrequency on the 
number of unknowns N of the FEM problem. 



wavelength of Ao = 489 nm is obtained for a nanowire length of L = 374 nm with a corresponding Q-factor of 
Q ~ 2.2. This is in very good agreement with the results of Section 13.31 Cross-sections through the 3D field 
distribution visualizing the resonant electric field distribution at L = 374 nm are shown in Figure [TOl (left). We 
have performed a convergence study for this specific cavity length. Here we investigate the convergence of the 
real and of the imaginary part of the complex eigenfrequency. For fixed geometrical parameters and for a finite 
element polynomial degree of p = 2 we increase the spatial resolution of the FEM mesh, leading to finer and finer 
meshes and FEM problems with increasing number of unknowns of the algebraic problem, N. Figure fTUI shows 
that for this specific setting, both, Ao and Q can be computed to a high numerical accuracy with a relative error 
below 10- 4 . 

4. CONCLUSION 

We have presented a numerical method for investigating nanoplasmonic devices like waveguides and cavities. 
Waveguide mode computations, 3D light scattering computations and 3D resonance mode computations have 
been performed. Accurate results have been obtained using higher-order finite elements and adaptive grid 
refinement. 
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